Correlations, fluctuations and stability of a finite-size network of coupled oscillators 
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The incoherent state of the Kuramoto model of coupled oscillators exhibits marginal modes in 
mean field theory. We demonstrate that corrections due to finite size effects render these modes 
stable in the subcritical case, i.e. when the population is not synchronous. This demonstration is 
facilitated by the construction of a non-equilibrium statistical field theoretic formulation of a generic 
model of coupled oscillators. This theory is consistent with previous results. In the all-to-all case, 
the fluctuations in this theory are due completely to finite size corrections, which can be calculated 
in an expansion in 1/N, where TV is the number of oscillators. The TV — > oo limit of this theory is 
what is traditionally called mean field theory for the Kuramoto model. 



I. INTRODUCTION 

Systems of coupled oscillators have been used to de- 
scribe the dynamics of an extraordinary range of phe- 
nomena including networks of neurons 0, 0], syn- 
chronization of blinking fireflies [J, [5| , chorusing of chirp- 
ing crickets @, neutrino flavor oscillations 0K arrays of 
lasers Q, and coupled Josephson junctions [3]. A com- 
mon model of coupled oscillators is the Kuramoto model 
[nil ], which describes the evolution of TV coupled oscilla- 
tors. A generalized form is given by 



(i) 



j 



where i labels the oscillators, 9i is the phase of oscillator 
i, f(9) is the phase dependent coupling, and the intrin- 
sic driving frequencies are distributed according to 
some distribution g{to). In the original Kuramoto model, 
f{8) = sin(0). Here, we consider / to be any smooth 
odd function. The system can be characterized by the 
complex order parameter 



'<*> = r(i)e^ (t) 



(2) 



where the magnitude r gives a measure of synchrony in 
the system. 

In the limit of an infinite oscillator system, Kuramoto 
showed that there is a bifurcation or continuous phase 
transition as the coupling K is increased beyond some 
critical value, K c ■ Below the critical point the steady 
state solution has r = (the "incoherent" state). Be- 
yond the critical point, a new steady state solution with 
r > emerges. Strogatz and Mirollo analyzed the lin- 
ear stability of the incoherent state of this system us- 
ing a Fokker-Planck formalism In the absence of 
external noise, the system displays marginal modes as- 
sociated with the driving frequencies of the oscillators. 
However, numerical simulations of the Kuramoto model 
for a large but finite number of oscillators show that the 
oscillators quickly settle into the incoherent state below 
the critical point. The paradox of why the marginally 
stable incoherent state seemed to be an attractor in sim- 
ulations was partially resolved by Strogatz, Mirollo and 



Matthews [12| who demonstrated (within the context of 
the TV — ► oo limit) that there was a dephasing effect akin 
to Landau damping in plasma physics which brought r to 
zero with a time constant that is inversely proportional to 
the width of the frequency distribution. Recently, Stro- 
gatz and Mirollo have shown that the fully locked state 
r = 1 is stable [HI but the partially locked state is again 
marginally stable Although dephasing can explain 
how the order parameter can go to zero, the question of 
whether the incoherent state is truly stable for a finite 
number of oscillators remains unknown. Even with de- 
phasing, in the infinite oscillator limit the system still has 
an infinite memory of the initial state so there may be 
classes of initial conditions for which the order parameter 
or the density exhibits oscillations. 

The applicability of the results for the infinite size 
Kuramoto model to a finite-size network of oscillators 
is largely unknown. The intractability of the finite size 
case suggests a statistical approach to understanding the 
dynamics. Accordingly, the infinite oscillator theories 
should be the limits of some averaging process for a finite 
system. While the behavior of a finite system is expected 
to converge to the "infinite" oscillator behavior, for a fi- 
nite number of oscillators the dynamics of the system will 
exhibit fluctuations. For example, Daido [H, [HI consid- 
ered his analytical treatments of the Kuramoto model 
using time averages and he was able to compute an ana- 
lytical estimate of the variance. In contrast, we will pur- 
sue ensemble averages over oscillator phases and driving 
frequencies. As the Kuramoto dynamics are determin- 
istic, this is equivalent to an average over initial phases 
and driving frequencies. Furthermore, the averaging pro- 
cess imparts a distinction between the order parameter 
Z and its magnitude r. Namely, do we consider (Z) or 
(r) = (\Z\) to be the order parameter? This is important 
as the two are not equal. In keeping with the density as 
the proper degree of freedom for the system (as in the in- 
finite oscillator theories mentioned above) , we assert that 
(Z) is the natural order parameter, as it is obtained via 
a linear transformation applied to the density. 

Recently, Hildebrand et al. [13] produced a kinetic 
theory inspired by plasma physics to describe the fluctu- 
ations within the system. They produced a Bogoliubov- 
Born-Green-Kirkwood-Yvon (BBGKY) moment hicrar- 
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chy and truncation at second order in the hierarchy 
yielded analytical results for the two point correlation 
function from which the fluctuations in the order param- 
eter could be computed. At this order, the system still 
manifested marginal modes. Going beyond second or- 
der was impractical within the kinetic theory formalism. 
Thus, it remained an open question as to whether going 
to higher order would show that finite size fluctuations 
could stabilize the marginal modes. 

Here, we introduce a statistical field theory approach 
to calculate the moments of the distribution function gov- 
erning the Kuramoto model. The formalism is equivalent 
to the Doi-Peliti path integral method used to derive sta- 
tistical field theories for Markov processes, even though 
our model is fully deterministic [H, ES EO, EH- The 
field theoretic action we derive produces exactly the same 
BBGKY hierarchy of the kinetic theory approach [It) ]. 
The advantages of the field theory approach are that 1) 
difficult calculations are easily visualized and performed 
through Fcynman graph analysis, 2) the theory is eas- 
ily extendable and generalizable (e.g. to local coupling), 
3) the field theoretic formalism permits the use of the 
renormalization group (which will be necessary near the 
critical point), and 4) in the case of the all-to-all ho- 
mogeneous coupling of the Kuramoto model proper, the 
formalism results in an expansion in 1/N and verifies 
that mean field theory is exact in the N — > oo limit. We 
will demonstrate that this theory predicts that finite size 
corrections will stabilize the marginal modes of the infi- 
nite oscillator theory. Readers unfamiliar with the tools 
of field theory are directed to one of the standard texts 
[22T | . A review of field theory for non-equilibrium dynam- 
ics is [13. 

In section [Til we present the derivation of the the- 
ory and elaborate on this theory's relationship to the 
BBGKY hierarchy. In section [TTTl we describe the com- 
putation of correlation functions in this theory and, in 
particular, describe the tree level linear response. This 
will connect the present work directly with what was 
computed using the kinetic theory approach [I?) • In sec- 
tion [iVl after describing two example perturbations, we 
calculate the one loop correction to the linear response 
and demonstrate that the modes which are marginal at 
mean field level are rendered stable by finite size effects. 
In addition, we demonstrate how generalized Landau 
damping arises quite naturally within our formalism. We 
compare these results to simulations in section IVl 



II. FIELD THEORY FOR THE KURAMOTO 
MODEL 

The Kuramoto model ((T|) can be described in terms of 
a density of oscillators in 9, u> space 



that obeys the continuity equation 



drj di] 
~di +UJ d6 
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f(6' - u/, t)rj(6, w, t)d9'duj' 

= (4) 

Equation (j4|) remains an exact description of the dynam- 
ics of the Kuramoto model {1} [17j . Although equa- 
tion (0} has the same form as that used by Strogatz 
and Mirollo (24|, it is fundamentally different because 
solutions need not be smooth. Rather, the solutions of 
Eq. (fj| are treated in the sense of distributions as defined 
by Eq. ([3]). As we will show, imposing smooth solutions 
is equivalent to mean field theory (the infinite oscilla- 
tor limit). Drawing an analogy to the kinetic theory of 
plasmas, Eq. @ is equivalent to the Klimontovich equa- 
tion while the mean field equation used by Strogatz and 
Mirollo [24| is equivalent to the Vlasov equation [25|, [2(| . 

Our goal is to construct a field theory to calculate the 
response functions and moments of the density rj. Even- 
tually we will construct a theory akin to a Doi-Peliti field 
theory [l8|, H3, HH , a standard approach for reaction- 
diffusion systems. We will arrive at this through a con- 
struction using a Martin-Siggia-Rose response field [27j . 
Since the model is deterministic, the time evolution of 
r]{9, u>, t) serves to map the initial distribution forward in 
time. We can therefore represent the functional proba- 
bility measure V[rj(9,ui,t)] for the density rj(9,ui,t) as a 
delta functional which enforces the deterministic evolu- 
tion from equation Q along with an expectation taken 
over the distribution Pofao] of the initial configuration 
r]o(9, m) = r](9,uj,to). We emphasize that no external 
noise is added to our system. Any statistical uncertainty 
completely derives from the distribution of the initial 
state. Hence we arrive at the following path integral, 

V[r)(6,uj,t)] 

■DrioV [m]S [N {C(r)(6, u, t)) - S(t - t )ri (6, lo)}} 

(5) 

The definition of the delta functional contains an ar- 
bitrary scaling factor, which we have taken to be N. 
We will show later that this choice is necessary for the 
field theory to correctly describe the statistics of 77. The 
probability measure obeys the normalization condition 

i = jv v vm : 

We first write the generalized Fourier decomposition of 
the delta functional. 

P[r,(9,u,t)} 



1 N 



T>fjVr)oV [rio] 



x exp -N / d6dudtrj[C(ri) - 8{t - t )r] {6,uj\ 



(6) 
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where fj(6, u>, t) is usually called the "response field" after 
Martin-Siggia-Rose [27j and the integral is taken along 
the imaginary axis. It is more convenient to work with 
the generalized probability including the response field: 

^[77,77] = J ZtyoPoMexp (-N J d9dudtrjC(ri) 



+ N / d9dLuij(9,uj,t )r] (9,uj) 



which obeys the normalization 



f 



Vr)T>fjV[r],fj;r)o] 



(7) 



(8) 



We now compute the integral over 770 in equation ijTj). 
which is an ensemble average over initial phases and driv- 
ing frequencies. We assume that the initial phases and 
driving frequencies for each of the N oscillators are inde- 
pendent and obey the distribution po(9,u>). V[r]o} repre- 
sents the functional probability distribution for the ini- 
tial number density of oscillators. Noting that rj (6, to) is 
given by 



770(0, w) = -^J2> 



(9) 



and J T>t]qVq[i']o\ = JY[id9idu!iPo(6i,u>i), one can show 
that the distribution from equation ([7]) is given by 



po(9,u) 



P [77,77] = exp [-N / d6dujdtf]C{ri) 



+ N\n{l + d9duj 



(10) 



In deriving Eq. (TIT)]) , we have used the fact that 



f^o^ofro] exp N / d6dujrj(0,u;,to)r) O (9,uj' 



J Y\ d9 i duj l p (9 i , uji) cxp ( ^ fj(6i,Wi,t ) 



N 



d9dLupo{6,Luy 1{8 > w > to) j 

= exp (n hi 1 + J d9dwp (9, w) (e^ e ' u ' to) - l) ^ 

(11) 

We see that the fluctuations (i.e. terms non-linear in 77) 
appear only in the initial condition of (|10[) . which is to 
be expected since the Kuramoto system is deterministic. 
In this form the continuity equation ([¥]) appears as a 
Langevin equation sourced by the noise from the initial 
state. 

Although the noise is contained entirely within the ini- 
tial conditions, it is still relatively complicated. We can 



simplify the structure of the noise in (|10p by performing 
the following transformation plj : 



(p(9,w,t) = nexp(-fj) 
(p(9,u>,t) + l = cxp(?7) 



(12) 



Under the transformation ([12]) , V\n,fi\ becomes V[<p,(p\, 
which is given by 



P[<p,<p] = exp(-NS[tp,tp\) 
where the action S[ip, 0\ is 



(13) 



S[<p, 0] = dtod9dt 











- K / dJd& {(p'tp + <p) — {f(9' - 6)<p'<p} 



111 



1 



d9dutp(9, u>, to)po(9, uS) 



(14) 



The form (TTJJ is obtained from the transformation (fTJ 
only after several integrations by parts. In most cases, 
these integrations do not yield boundary terms because 
of the periodic domain (i.e. those in 9). In the case of the 
dt operator, however, we are left with boundary terms of 
the form [ln(</5 + 1) — l]<ptp. These terms will not affect 
computations of the moments because of the causality of 
the propagator (see section [Til A[) . 

We are interested in fluctuations around a smooth so- 
lution p(6,u),t) of the continuity equation ((4]) with initial 
condition p{6, to, to) = po(9,u). We transform the field 
variables via ^ — ip — p and ip = (p in (|14|) and obtain the 
following action: 



Sty, 0] 



dtudOdtip 



d d , 

di + L °d9^ 



+ kJ du'dS'^ {f (9' - 9) W P + p'tjj + W)} 
+ K [ Aoj'd&i,'^ {f(9' - 9) W + p') (yj + p)} 



E 

k=2 



(_l)fe+l 



d9duj r i/j(9, to, to)po(9, oj) 



(15) 



For fluctuations about the incoherent state: p(9,u,t) = 
Po(6,oj) = g(u))/2ir, where g(u>) is a fixed frequency dis- 
tribution. The incoherent state is an exact solution of 
the continuity equation ([4]). Due to the homogeneity in 
9 and the derivative couplings, there are no corrections to 
it at any order in 1/N. The action (fTS"]) with p = g(ui)/2ir 
therefore describes fluctuations about the true mean dis- 
tribution of the theory, i.e. (r)(9,ui,t)) = g(uj)/2w. We 
can evaluate the moments of the probability distribution 
(fT3"| with (fTBj) using the method of steepest descents, 
which treats 1/N as an expansion parameter. This is a 
standard method in field theory which produces the loop 
expansion [22]. We first separate the action (fl~5|) into 
"free" and "interacting" terms. 



Sty,^\ = S F ty^]+S I ty^} 



(16) 
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where 



x m 4 • x 



Po{x,t\x',t') 



S F [ip,ip] 

+ K J du'de'^- e {f{6> -e)W P + P >i>)} 



dxdtdx dt ip(x , t')ro(x, t; a; , t)ip(x, t) 



(17) 



cLodOdtip 



du'de'-^me'-e) (W)} 



dw'de'i*'-^ {f(0' - e) w + p') a> + P )} 



-E 



fe+i r 



d9dtuip(9, u), to)po(9, 



(18) 



In deriving the loop expansion, the action is expanded 
around a saddle point, resulting in an asymptotic series 
whose terms consist of moments of the Gaussian func- 
tional defined by the terms in the action (|15p which are 
bilinear in ip and ip, i.e. V']- Hence, the loop expan- 

sion terms consist of various combinations of the inverse 
of the operator r , defined by Sp, called the bare prop- 
agator, with the higher order terms in the action, called 
the vertices. Vertices are given by the terms in Si. 

The terms in the loop expansion are conveniently rep- 
resented by diagrams. The bare propagator is repre- 
sented diagrammatically by a line, and should be com- 
pared to the variance of a gaussian distribution. Each 
term in the action (other than the bilinear term) with n 
powers of ip and m powers of ip is represented by a vertex 
with n incoming lines and m outgoing lines. The initial 
state vertices produce only outgoing lines and, like the 
non-initial state or "bulk" vertices, are integrated over 
6, u), and t for each point at which the operators are de- 
fined. The bulk vertices are represented by a solid black 
dot (or square, see Figure [T]) and initial state vertices 
by an open circle. The bare propagator and vertices are 
shown in Figure [T] Unlike conventional Fcynman dia- 
grams used in field theory, the vertices in Figure [1] rep- 
resent nonlocal operators defined at multiple points. In 
particular, the initial state terms involve operators at a 
different point for each outgoing line. Although uncon- 
ventional, this is the natural way of characterizing the 
1/N expansion. 

Adopting the shorthand notation of x = {9,cu}, each 
arm of a vertex must be connected to a line (propagator) 
at either x or x and lines connect outgoing arms in one 
vertex to incoming arms in another. The moment with 
n powers of ip and m powers of ip is calculated by sum- 
ming all diagrams with n outgoing lines and m incoming 
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FIG. 1: Diagrammatic (Feynman) rules for the fluctuations 
about the mean. Time moves from right to left, as indicated 
by the arrow. The bare propagator Po(x, t\x' , t') (see Eq. (|31[) ) 
connects points at x' to x, where x = {8, ui}. Each branch of 
a vertex is labeled by x and x' and is connected to a factor of 
the propagator at x or x' . Each vertex represents an operator 
given to the right of that vertex. The ". . . " on which the 
derivatives act only include the incoming propagators, but 
not the outgoing ones. There are integrations over 9,9',u),u)' 
and t at each vertex. 



lines. This means that each diagram will stand for sev- 
eral terms which are equivalent by permutations of the 
vertices or the edges in the graph, equivalently permuta- 
tions of the factors of ip and ip in the terms in the series 
expansion. In a typical field theory, this results in combi- 
natoric factors. In the present case, diagrams which are 
not topologically distinct can produce different contribu- 
tions to a given moment. Nonetheless, we will designate 
the sum of these terms with a single graph. Generically, 
the combinatoric factors we expect are due to the ex- 
change of equivalent vertices, which typically cancel the 
factorial in the series expansion. Additionally, each line 
in a diagram contributes a factor of 1 /N and each vertex 
contributes a factor of N. Hence each loop in a diagram 
carries a factor of 1/N. The terms in the expansion with- 
out loops are called "tree level". The bare propagator 
is the tree level expansion of (ip(9,id,t)ip(9',w',t')). The 
tree level expansion of each moment beyond the first car- 
ries an additional factor of 1/N, i.e. the propagator and 
two-point correlators are each 0(1/N). 

Mean field theory is defined as the N oo limit of this 
field theory. In the infinite size limit, all moments higher 
than the first are zero (provided the terms in the series 
are not at a singular point, i.e. the onset of synchrony). 
Hence, the only surviving terms in the action (|15|) are 
those which contribute to the mean of the field at tree 
level. These terms sum to give solutions to the continuity 
equation (|4|). If the initial conditions are smooth, then 
mean field theory is given by the relevant smooth solution 
of (j4]). In most of the previous work (e.g. [l2|, H3]), 
smooth solutions to ^ were taken as the starting point 
and hence automatically assumed mean field theory. 

We can now validate our choice of TV as the correct 
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scaling factor in the delta functional of ([5]) by considering 
the equal time two-point correlator. Using the definition 
of r\ from ((3]) we get 

(r](x,t)ri(x',t)) = C(x,x',t) + p(x,t)p(x',t) 

- Trp(x,t)p(x',t) + —S(x-x')p(x',t) 

(19) 

where p(x,t) = (r)(x,t)) = (<p(x,t)). Using the fields ip 
and ip (defined in (|12[)) and taking r\ at different times 
gives 

(■n(x,t)ri(x',t')) 
= ([tp(x, t)ip(x, t) + <p(x, t)} [<p(x>, t'Mx',t') + cp(x', t')\) 

(20) 

for t > t'. The response field has the property that ex- 
pectation values containing (p(x, t) are zero unless an- 
other field insertion of <p{x,t) is also present but at a 
later time (this is because the propagator is causal; it is 
zero for t — t' < 0). Therefore 

(r](x,t)r){x\t')) = (p(x,t)p(x' ',*')) 

+ (p(x,t)p(x',t'))(p(x',t')) (21) 

As we will show later when we discuss the propagator in 
more detail, we have 

Km(<p(x,t)<p(x'J)) = ±6(x-x') (22) 

Comparing (|2"Tj) in the limit t — > t' with (|19[) allows the 
immediate identification of 

( V {x,t)ip{x' ,t)) = C(x,x',t) + p(x,t)p(x',t) 

- ^p(x,t)p(x',t) (23) 

C(x,x',t) is the two oscillator "connected" correlation 
or moment function. This is consistent with (|15|) which 
gives 

{ip(x,t )<p(x',t )) = p(x,t )p(x',t ) - ^-p(x,t )p(x',t ) 

(24) 

as the initial condition. Thus, comparing the second mo- 
ment of rj using the Doi-Peliti fields (f2"0|) with the expres- 
sion from the direct computation given by Q19[) shows 
that the factor of N in the delta functional of ([5]) was 
necessary to obtain the correct scaling for the moments. 

The Doi-Peliti action (fTS"|) can also be derived by con- 
sidering an effective Markov process on a circular lattice 
representing the angle 8 where the probability of an oscil- 
lator moving to a new point on the lattice is determined 
by its native driving frequency u> and the relative phases 
of the other oscillators (see Appendix [C]). This is the 
primary reason we refer to the theory as being of Doi- 
Peliti type. The continuum limit of this process yields a 



theory described by the action (fl~5|) . The Markov picture 
provides an intuitive description and underscores the fun- 
damental idea that we have produced a statistical theory 
obeyed by a deterministic process. 

Although our formalism is statistical, we emphasize 
that no approximations have been introduced. The sta- 
tistical uncertainty is inherited from averaging over the 
initial phases and driving frequencies. This formalism 
could be applied to a wide variety of deterministic dy- 
namical systems that can be represented by a distribu- 
tional continuity equation like Eq. (f5|) . In general, a solu- 
tion for the moment generating functional for our action 
(fT5")) is as difficult to obtain as solving the original sys- 
tem. The advantage of formulating the system as a field 
theory is that a controlled perturbation expansion with 
the inverse system size as the small parameter is possible. 



A. Relation to Kinetic Theory and Moment 
Hierarchies 

The theory defined by the action (|15p is equivalcntly 
expressed as a Born-Bogoliubov-Green-Kirkwood-Yvon 
(BBGKY) moment hierarchy starting with the continuity 
equation ((4]). To construct the moment hierarchy, one 
takes expectation values of the continuity equation with 
products of the density, tj(8, uj, t). This results in coupled 
equations of motion for the various moments of r/(8, u, t), 
where each equation depends upon one higher moment in 
the hierarchy. The Dyson-Schwinger equations derivable 
from (fT5|) with (p) = are exactly the BBGKY hierarchy 
derived in Ref. [17|. Thus the kinetic theory and field 
theory approaches are entirely equivalent. 

The moments of rj can be computed in the BBGKY hi- 
erarchy by truncating at some level. In Ref. [T3|, this was 
done at Gaussian order by assuming that the connected 
three-point function was zero. The first two equations 
of the hierarchy then form a closed system which can by 
solved. The first two equations of the BBGKY hierar- 
chy [13] are 

dp dp 
~dt +UJ d8 

+ K—J f(9'-8)p(x,t)p(x',t)d6'dLj' 

= - K ggj J q f(0'-8)C(x,x',t)d6'dLj' (25) 

where p(x,t) = (j](x,t)}, and the connected (equal-time) 
correlation function 

C(x,x',t) = (rj(x,t)r](x',t)) - p(x,t)p(x' ,t) 

+ jjPfa t)p(x',t) - —S(x - x')p(x', t) 

(26) 

obeys 
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/oo pin a poo pin a 

/ ^ r f(03-0 1 )p(x 1 ,t)C(x 2 ,x 3 ,t)d9 3 diJ 3 + K / f(e 3 -e 2 )p(x 2 ,t)C(x 3 ,x u t)de 3 cL; 3 } 

= - + - e 2 )]p(x 1 ,t) P {x 2 ,t). 

(27) 



In the field theoretic approach, instead of truncating 
the BBGKY hierarchy, one instead truncates the loop 
expansion. Truncating the moment hierarchy at the mth 
order is equivalent to truncating the loop expansion for 
the ith moment at the (m — i)th order. Thus the solu- 
tion to the moment equations ([25]) and ([27]) is the one 
loop expression for the first moment and the tree level 
expression for the second moment. The advantage of us- 
ing the action (|15p is that the terms in the perturbation 
expansion are given automatically by the relevant dia- 
grams at any level of the hierarchy. Ref. [TtJ suggested 
that a higher order in the hierarchy would be necessary to 
check whether the mean field marginal modes are stable 
for finite N. We demonstrate below that the field the- 
ory facilitates the calculation of the linearization of equa- 
tion (|23|) to higher order in 1 /N and show that marginal 
modes are stabilized by finite size fluctuations. 

One can compare this approach to the maximum en- 
tropy approach of Rangan and Cai [28[ for developing 
consistent moment closures for such hierarchies. In the 
moment hierarchy approach of Ref. [l7j | , moment closure 
is obtained via the somewhat ad hoc approach of setting 
the nth cumulant to zero. In contrast, Rangan and Cai 
maximize the entropy of the distribution subject to cer- 
tain normalization constraints. The moment closure is 
facilitated by constraining higher moments from the hi- 
erarchy. However, one still must solve the resulting equa- 
tions. In the loop expansion, moment closure is obtained 
implicitly via truncating the loop expansion. The loop 
expansion approach offers the advantage of providing a 
natural means for determining when the approximation, 
thus the implicit closure, breaks down and avoids deal- 
ing with the moment hierarchy explicitly. In fact, Ran- 
gan and Cai's procedure has a natural interpretation in 
field theory, namely the minimization of a generalized 
effective action in terms of various moments. The sim- 
plest and most common is the effective action in terms 
of the mean field, which is the generating functional of 
one partical irreducible (1PI) graphs [22j|. The next level 
of approximation is a generalized effective action (the 
"effective action for composite operators" [2{|) in terms 
of the mean and the two-point function (or functions), 
which is the generating functional of two particle irre- 
ducible (2PI) graphs. One can continue in this way. The 
equations of motion of these effective actions will produce 



a closure of the moment hierarchy implicit in the action 
for the theory. At tree level these equations will be equiv- 
alent to those produced by Rangan and Cai's maximum 
entropy approach. The loop expansion allows for sys- 
tematic corrections to these equations without explicitly 
invoking higher equations in the hierarchy. 



III. TREE LEVEL LINEAR RESPONSE, 
CORRELATIONS AND FLUCTUATIONS 

As a first example, we reproduce the calculation of the 
variation of the order parameter Z , which was calculated 
previously using the BBGKY moment hierarchy [T3] ■ To 
do so requires the calculation of the tree level linear re- 
sponse or bare propagator and the tree level connected 
two-oscillator correlation function. 



A. The Propagator 

The propagator P(0, w, t\9' , w', t') is given by the ex- 
pectation value 

P(0,(j,t\e',L>',t') = (<p{6,uj,t)(p{6',u)\t')) (28) 

It is the linear response of (@|. This can be shown by 
considering a small perturbation Spo to the initial state 
p in the action (| 1 5|) . Expanding to first order then yields: 

Sp(6,Lj,t)=N J dO'duj' (tp(6,Lj,t)<£(6',u',t'))5p (6',ij') 

(29) 

The tree level linear response or bare propagator 
Pq(9, u>, t\6' , ui', t') = Pq(x, t; x', t') is the functional in- 
verse of the operator Tq defined by the free part of the 
action (|17[) . The bare propagator is therefore given by 

M 

T -P ee J dx"dt"T (x,t;x",t")P (x",t";x',t') 

= ±5(x-x')5(t-t') (30) 
Using the action (|15|) with Eq. (|3"0|) gives 
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+ K—J / f(0 1 -9) P (x,t)P o (x 1 ,x',t-t')de 1 du 1 



P (x,x',t-t') 



—8{9-8')8{to-uj')5(t-t') 



(31) 



Due to the rotational invariance in 9 of f(9), 
P(x, t; x', t') = P{8 - 9', to, to', t - if). 

In the incoherent state, p{9,uj,t) = g(uj)/2ir. Thus, for 
f{6) odd, Eq. dHTJ becomes 



K 



d d 

¥t +UJ dl 

2tt 961 



Po(£, ^ j £ ^ t ) 

OO />27T 



f(9 1 -9)P (x 1 ,x',t-t')de 1 dto 1 



1 



= ^(9-^f(w-i/){({-*') 



(32) 



We can invert this equation using Fourier and Laplace 
transforms. 

Taking the Fourier transform of Eq. (|3"2")) (with respect 
to 9 and 9') yields 



d_ 

at 



Po(n, to; m, to', t — t') 



+ inKg(to)f(—n) j Po(n, toy, m, u>', t — t')dui 
1 



2-kN 



S(t - t')S(to - uj')8 n+m 



(33) 



where we use the following convention for the Fourier 
transform 



i 

2^ 



2tt 



f{oy 



f(9) = £/(n)e 



illO 



(34) 



Hereon, we will suppress the index m since the propaga- 
tor must be diagonal (i.e. Po(n, m)) oc 5 m + n ). 
We Laplace transform in r = t — t' to get 

[s + mcj] Po(n, a/, s) 

/OO 
P (n, uji, uj', s)du 1 
-OO 

1 



2itN 

using the convention 



6(u) — to') 



(35) 



f{t)e~ ST dT 



1 



2,771 ' 



(36) 



where the contour C is to the right of all poles in f{s). 

We can solve for Po(n, u, to', s) using a self-consistency 
condition. Integrate (|35[) over to after dividing by s + mcj 
to get 



duoPoin, uj, uj' , s) 

inKg(u))f{-n) , 

duj / duJiPo(n,uJi,uj , s) 

s + inuj 

1 1 

27T./V s + inuj' 

which we can solve to obtain 



(37) 



. . 1 1 1 

dujP [n,uj,uj ,s) = — -— — (38) 

2nN s + itvjj' A n (s) 



where 



A„(s) = 1 + inKf(-n) / du> 



s + into 



(39) 



A n (s) is defined for Re(s) < via analytic continuation. 
In the kinetic theory context of an oscillator density obey- 
ing the continuity equatio n ffl , A„(s) is analogous to a 
plasma dielectric function |l7| . If we assume that g(to) is 
even and f(9) is odd, then there is a single real number 
s n such that A n (s n ) = 0. (Mirollo and Strogatz proved 



that there is at most one single, real root of (|39|) and that 
it must satisfy Re(s) > [0,|3fl|. In our case, A n (s) is 
defined for Re(s) < not by (f3"9"| . but rather via ana- 
lytic continuation.) Using (|39|) in (|35[) and solving for 
P (n,uj,to',s) gives 



P Q (n,uj,uj',s) 



1 S(uj-uj') 
2nN s + into 

1 inKg(uj)f{-n) 1 
2k N (s + into) (s + into') A n (s) 



(40) 

Here we identify the spectrum with the zeroes of To or, 
cquivalently, the poles of the propagator, as these will de- 
termine the time evolution of perturbations. Analogous 
to the analysis of Strogatz and Mirollo [24[ we define the 
operator O by 

O[b n (to,t)] = inub„(co,t) 

/>oo 

+ inKf- n I b n (uj 1 ,t)g(to 1 )dto 1 (41) 
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a) 

> 



b) 




FIG. 2: Diagrams for the connected two-point function at tree 
level a) and to one loop b). 



This is equivalent to ([26)) . which was computed in 
Ref. [ljj when t\ = t 2 . If the initial phases are un- 
correlated (i.e. C(x±, 0; x 2 , 0) = 0) then at tree level 
C{x\, ti; X2, t 2 ) is given by the diagram shown in Fig- 
ure [2] a). It is comprised of vertex III (see Figure [1} com- 
bined with a bare propagator on each arm. For general 
(i.e. not odd) f(8), the diagram in Figure [5] a) actually 
corresponds to two different terms because the arms of 
vertex III can be interchanged, giving two terms in equa- 
tion l|43p below. Unlike, conventional field theory these 
interchanges are not symmetric. These two terms are 
equal when f(9) is odd. More generally other vertices 
do not exhibit the symmetries typical of Feynman dia- 
grams even for odd f(0). Applying the Feynman rules 
then gives at tree level 



(cf. equation (|62[) below). The continuous spectrum of 
O consists of the frequencies into whereas the discrete 
spectrum (according to Ref. [24[) only exists for K > 
K c . Consistent with that approach (i.e. linear operator 
theory) . we identify the poles in P due to s + inu> as the 
continuous spectrum and those due to the zeroes of A n (s) 
as the discrete spectrum. If A„(s) is not analytically 
continued for Re(s) < 0, it will not have zeroes for that 
domain, as in Ref. [24|. However, zeros can exist for 
Re(s) < when A„(s) is analytically continued and this 
is why analytic continuation is of such crucial importance 
to the conclusions of Ref. [l2| . 



C(xi,t 1 ;x 2 ,t 2 ) 
x P (x 1 ,x' 1 ,t 1 - t')P (x 2 ,x' 2 ,t 2 - t') 

x i^m 0i) + JLm 2 )]<7K)<?k) (43) 



B. Correlation Function 

The connected correlation function (cumulant func- 
tion) is given by 

C(xi,t 1 ;x 2 ,t 2 ) = (^(x 1 ,t 1 )i>(x 2 ,t 2 )) (42) 



where t = min(ti, t 2 ). This is essentially identical to the 
ansatz used in Ref. [17| for the solution of the second 
moment equation in the BBGKY hierarchy (|2"T|) . For 
f(9) = sin# and t\ > t 2 , Fourier transforming Eq. (|4"5|) 
and inserting the bare propagator from Eq. (|4H|) (after 
an inverse Laplace transformation) gives 



K 1 



g(u\)g{u2) 



(i Wl + ^-f)(-i W2 + ^-f) 



,(i(u>i—u>2)t 2 ) 



1 



- lu 2 ) i{w\ - uj 2 ) 
K (ilu 2 - 3») 



3 iwi(ti-* 2 ) 



2 (K_^_ ilVl){{ ^_K ) 2 + {uj2) 2 ) 
K (-tW! - ^f) 



^(^-f+i^2)(t 2 ) _ g-(^-f )(tl- 



t 2 ) 



2 (f _^ +iw2 )(( 



#) 2 + M 2 ) 



-(-f—-2-iU>l)t 2 _ I 



-iu 1 (t 1 —t 2 ) 
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K 2 



A(K - 2^) (f - ^ - i Wl )(f - 4f + iw a ) 



(e-^-^* 2 - l) 



,-(K c -Jf))(ti-i a ) 



The equal time correlator is given by (t\ = t 2 = r): 



(44) 



A' 1 

Ci(wi,w 2 ,t) = — -^g(wi)g(w 2 ) 



(iW! + ^)(-^ 2 + ^) 



(icJi + ^ - §)(-ioj 2 + %f - f ) V - w 2) - w 2 ) 



(^2 " ^) 



2 (f -^- iWl )(( i |--f) 2 + (^) 2 ) 
£ HtJi - ^ ) 

2 (f _^+^ 2 )((^_f)2 + (a;i) 2 ) 
-2 



2 2 

1 



,-(-2^- — +iw 2 )r 



-(^-f-^i)^ _ 1 



4(K - 2f) (f - ^ - i Wl )(f - ^ + i«a) 



-(_ftT c -K)r 



- 1 



(45) 



C_i = C* and the other modes vanish. Note that the 
initial condition C\(lo,lo' , 0) = is satisfied and that the 
time constants and frequencies which appear are every 
possible way of pairing those from the tree level propa- 
gator. 

For illustrative purposes. Figure [2] b) shows the one 
loop diagrams which contribute to C; these diagrams are 
0(1/N 2 ). We should note here the special role played by 
the diagrams with initial state terms, in particular the 
vertex proportional to ip(0, u>, to) 2 - This diagram evalu- 
ates to exactly the same result as (|4"4"]l . with an additional 
factor of —1/N. It serves to provide the proper normal- 
ization for the two point function, which should go as 
(N — l)/N since the self-interaction (diagonal) terms are 
not included. The other diagrams diverge faster as one 
approaches criticality (K = K c ). They are of negligible 
importance at small coupling but become increasingly 
important near the onset of synchrony. 



C. Order Parameter Fluctuations 

We now compute the fluctuations in the order parame- 
ter Z given in Eq. ^ . The variance of Z (second moment 
(ZZ)) is given by 



(ZZ) = (r 2 (t)) 



= / dwduj'd0d0> '{< n {uj,6,t)<n{J ,6' ',*)) 



Using equation (flU|) in equation (|4"6"| gives 



' ft' t\\A - s ') 



(46) 



{r 2 (t)) = j dujdLu , d0d0'C(x 1 x\t)e l{e ' e ' ) + i (47) 



since in the incoherent state p(x,t) = g(ci>)/2ir is inde- 
pendent of 8, so that (Z) = 0. Hence 



(r 2 (r)) = 4tt 2 / dwdu'C-t^, u', r) + 



N 



(48) 



which evaluates to ITi 



! (r)> 



iKNir 
x Res 



da 



AiOO 



Ai(f - go) - 1 
Ai(s - s ) 

1 

N 



-e 

s 



77 (49) 



The time evolution of (r 2 ) is then determined by the 
poles of Ai(s). As an example, consider f(9) = s'm9 and 
g(us) a Lorentz distribution (see Appendix |B|) : we have 
the result 



(r 2 (r)) 



1 K r 



1 K 



(K c -K) 



N K r -K N K r - K 



(50) 



where K c = 27. Note that this diverges as r — > 00 for 
K = K c . In the mean field limit N — > 00, (?' 2 ) = as 
expected. As was shown in Ref [lT}, the tree level calcu- 
lation adequately captures the fluctuations except near 
the onset of synchrony (K = K c ). An advantage of the 
field theoretic formalism is that it allows us to approach 
even higher moments without needing to worry about 
the moment hierarchy. In particular, for f(8) = sm9 it 
is straightforward to show that higher cumulants, such 
as ((ZZ) 2 } — (ZZ) 2 , must be zero at tree level because 
of rotational invariance (more precisely, any cumulant of 
Z higher than quadratic). The cumulants are given by 
graphs which are connected. Vertex III produces two 
lines with wave numbers n = ±1. Additionally, the IV 
and V vertices impose a shift in wave number, whereas 
Z and Z project onto ±1. In order to calculate these 
higher fluctuations it is necessary to go to the one loop 
level. Note that this does not imply that the higher cu- 
mulants of i] are zero. Figure [5] b) gives the diagrams 
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for the correlation function at one loop. The one loop 
calculation would also give a better estimate for (r 2 ), es- 
pecially nearer to criticality. The non-interacting distri- 
bution is Gaussian, with non-Gaussian behavior growing 
as one approaches criticality. 



IV. LINEAR STABILITY AND MARGINAL 
MODES 

We analyze linear stability by convolving the linear 
response or propagator with an initial perturbation. Al- 
though we are free in our formalism to consider arbitrary 
perturbations, we will consider two specific kinds for il- 
lustrative purposes. In the first case we perturb only the 
angular distribution. In the second case we consider a 
perturbation which fixes one oscillator to be at a given 
angle and frequency a; at a given time t. We first calcu- 
late the results at tree level which reproduces the mean 
field theory results of Ref. [24| . In particular, we arrive 
at the same s pec trum and Landau damping results of 
Refs. [HI and [l2| . We then define and calculate the op- 
erator r (an extension of To) to one loop order which 
ultimately allows us to calculate the corrections to the 
spectrum to order 1/N. 



A. Mean field theory 

The bare propagator which is the full linear response 
for mean field theory is given by Eq. (|40|) . The zeroes 
of the operator T with respect to s specify the spectrum 
of the linear response. There is a set of marginal modes 
(continuous spectrum) along the imaginary axis spanning 
an interval given by the support of g{oj). There is also a 
set of discrete modes given by the zeros of the dielectric 
function A„(s) as was found in Rcf. [24j . aside from the 
issue of the analytic continuation of A„(s). 

However, even though there are marginally stable 
modes, the order parameter Z can still decay to zero due 
to a gen eralized Landau damping effect as was shown in 
Ref. [l2|. Consider a generalization of the order param- 
eter 



*«(*) = ^£ 



(51) 



which represents Fourier modes of the density integrated 
over all frequencies: 



Z n (t) = / d9duji](6,uj,t)e 



and hence 



(52) 



(Z n {t)) = J d9dujp(9,uj,t)e ine 

= 2tt / dujp(—n 7 ui 7 t). (53) 



In the incoherent state, p(0,u>,t) is independent of 9 so 
(Z n ) is zero for n > 0. The density response Sp(0,u),t) 
to an initial perturbation 5p(0,w,O) is given by 

/OO />2"7T 
/ P o (x,x',t)6p(0',uj',O)d0'du>' 
-oo JO 

(54) 

Recall from the definition of the action ([15]) . the prop- 
agator operates on an initial condition defined by Npo- 
The perturbed order parameter thus obeys 



5(Z n (t)) = / d0duj6p(0 1 uj 1 t) 



AnO 



(55) 



We will show that for any initial condition involving a 
smooth distribution in frequency and angle, 5(Z n (t)) will 
decay to zero. However, for non-smooth initial perturba- 
tions, 8(Z n (i)) will not decay to zero but will oscillate. 
We first consider an initial perturbation of the form 



where 



6p(0,Lo,O)=g(Lo)c(8) 



c{0)d0 = 



(56) 



(57) 



Inserting into (p3[) yields 

/ P o (x,x',t)c{0')g{w')d6'du' 

-oo JO 

(58) 

which is consistent with the perturbation considered in 
Ref. [24[. Taking the Laplace transform of (f5"8"|) gives 

/•OO 

5p n (uj, s) — Nc n / Po(n, w, u>', s)g{w')duj l (59) 



Using the tree level propagator (|4TJ|) , we can show that 

1 g(u) 1 



P (n,uj,uj',s)g(uj')duj' 

ZirlM s + into A n (s) 
where A„(s) is given in (|39|) . Hence 

c„ g(uj) 1 



5p(n,w,s) 



2n s + inuj A„(s) 



(60) 



(61) 



From Eq. (foTj) , we see that the continuous spectrum is 
given by into and the discrete spectrum by the zeros of 
A„(s). If we define b n (u>,t) = Sp(n,uj,t)/(Nc n g(uj)) then 
using (|58|) and (|40|) we can show 



dt 



+ inuj 



b„(cj,t) 



+ inKf(-n) j b n (uji,t)g(uji)dwi 
S(t) 



2nN 



(62) 
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which is equivalent to the linearized perturbation equa- 
tion derived by Strogatz and Mirollo [24| . with the ex- 
ception that ([62]) includes the effects of the initial config- 
uration through the source term proportional to S(t). 

Inserting (fBTj) into the the Laplace transform of (f53|) 
yields 



8(Z n (s)) = c_ n / 9 U du 



s — inuj A- n (s) 

A-n(a) ~ 1 l_ 

" -inKf(ri) A_„(s) ' 

For f(6) = sin (9, /(±1) = =r/ 2 which leads to 

A_i(«)-1 1 



S(Z 1 (a))=c- 



K/2 A_i(«)" 



(63) 



(64) 



We note that 8Zi(s) is identical to what was calculated 
in Ref. [l2j in which it was shown that 8Z\(t) ~ * as 
t — > oo. Even in the presence of marginal modes, the 
order parameter decays to zero through dephasing of the 
oscillators. This dephasing effect is similar to Landau 
damping in plasma physics. 

We can see this explicitly for the case of the Lorcntz 
distribution 



1 7 



7T 7 2 + LU 2 



From (|3"9"j) we can calculate 



A±i(s) = 



s + 7 



(65) 



(66) 



and A„(s) = 1 for n 7^ ±1. The zero of A±i(s) is at 
s±i = —(7 — K/2), which provides a critical coupling 



K c = 2 7 



(67) 



above which the system begins to synchronize. The in- 
coherent state is reached when K < K Cl which gives 
s±i < 0. Thus 



8(Z ±1 ) = c Tie -(^> 
8(Z n ^ ±1 ) = c_ n e-H^ 



(68) 



and so the initial perturbation is 



1 

N 



2tt 



■8(6 - e )5(uj - u ) 



(70) 



Inserting into (|54[) gives the time evolution of this initial 
perturbation 

8p(6, u,t) = -7^^ + P^, w, t; 0o, w 0l (71) 

Substituting into (|55|) and taking the Lapace transform 
gives 



8(Z n (s)) = 2tt J doj8p~- n (uj,s) 

= 2n / dLuP(—n,LU,u!o,s) 



N s — inojQ A_ n (s) 



(72) 



There are therefore two modes in 8Z n (t), one which de- 
cays due to dephasing (determined by the zero of A_„(s)) 
and one which oscillates at frequency luq. Thus, for this 
perturbation, the tree level prediction is that 8(Z) is not 
zero but oscillates. Inverse Laplace transforming l|72[). 
gives the time dependence of the order parameter 



8(Z 1 (t)) = 



1 



+ (7"f) 



7 - y J + - yiwo 



-icj + 7 



_ £ 1 ^ e -(7-4F )* 



(73) 



The perturbed oscillator has phase 9q — ^ot- It can al- 
ways be located; no information is lost as the time evo- 
lution progresses. Hence, for a single oscillator pertur- 
bation, the tree level calculation predicts that the order 
parameter will not decay to zero. In the next section, we 
show that to the next order in the loop expansion, which 
accounts for finite size effects, the marginal modes are 
moved off of the imaginary axis stabilizing the incoher- 
ent state, and the order parameter for a single oscillator 
perturbation decays to zero. 



Hence angular perturbations decay away in the order pa- 
rameter. 

Landau damping due to dephasing is sufficient to de- 
scribe the relaxation of Z n (t) to zero for a smooth per- 
turbation. However, for non-smooth perturbations, this 
may not be true. Consider the linear response to a stim- 
ulus consisting of perturbing a single oscillator to have 
initial position 8q and frequency too- 



(N -1)^ + 8(9 -e )8(oj- co) 
Ztt 



(69) 



B. Finite size effects 



Let us define a generalization of the operator Tq. 

r • P = ±8(x - x')8(t - t') (74) 

where P without a subscript denotes the full propagator 
and the operator V is the functional inverse of P. We 
can estimate the effect of finite size on the stability of the 
incoherent state by calculating the one loop correction to 
the operator T. We will see that the one loop correction 
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-<y<y 

FIG. 3: The diagrams contributing to the propagator at one 
loop order, organized by topology. We consider d) to be of 
different topology than c) because it is equivalent to a tree 
level diagram with an additional factor of 1/N, due to the 
initial state vertex. 

— ^ 

FIG. 4: Diagrammatic equation for the propagator. The dou- 
ble lines represent the summation of the entire series in 1/N 
for the propagator. 

produces the effect, among others, of adding a diffusion 
operator to T, which is enough to stabilize the continuum 
of marginal modes because the continuous spectrum is 
pushed off the imaginary axis by an amount proportional 
to the diffusion coefficient. 

We calculate the correction to T to one loop order. The 
propagator is represented by diagrams with one incom- 



ing line and one outgoing line. There are four groups of 
diagrams which contribute to the propagator at one loop 
order. They are shown in Figure [3] and are labeled by 
a), b), c), and d). Using these graphs to calculate the 
propagator to order 1/N is not sufficient to demonstrate 
the behavior of the spectrum to order 1/N. However, 
we can use these graphs to construct an approximation 
of r to order 1/N and derive the spectrum from this. 
If we denote the full propagator (i.e. the entire series 
in 1/N for P) by a double line, we can approximate the 
full propagator recursively by the diagrammatic equation 
shown in Figure 2] The only terms which are neglected 
in this relation are those which are from two loop and 
higher graphs and therefore would contribute 0(1/ N 2 ) 
to r. Readers familiar with field theory will note that we 
are simply calculating the two point proper vertex, which 
is the inverse of the full propagator, to one loop order. 

If we act on both sides of this equation with To (the 
operator whose inverse is the tree level propagator), we 
arrive at an equation of the form: 

r • P = ^S(x - x')S(t -i>)-Ti>P (75) 

where we have implicitly defined the one loop correction 
to r, which we label Ti. The action of To converts the 
leftmost propagator in each diagram into a delta func- 
tion, so that the delta function term in (|75p arises from 
the tree level propagator line and Ti is then comprised of 
the loop portions of the remaining diagrams (the "ampu- 
tated" graphs). We denote the contribution to Ti from 
each group of one loop diagrams by Ti r (9, u>; 4>,rj;t — t'), 
where r represents a, b, c, or d indicating the group of 
diagrams in question. Tu, = because the derivative 
coupling acts on the incoherent state, which is homoge- 
neous in 9. The equation of motion for the one loop 
propagator P\(x, x', t) then has the form 

r • Pr + r x ■ Pr = -^6(6 - 6')5(u - u')5(t - t 1 ) (76) 
where To ■ P\ is given by 



and 



P\(x, X , t — t) 



K 



d_ 

de 



oo JO 



f{6 x - 6)p(x,t)P 1 (x- L ,x',t- t^dOxdw-i 



(77) 



Ti • Pi = 



2tt roc 



drj / dt"T la (6,u ] cj), m t-t")P 1 ((l),r i ,t" ] 6 , ,u>'-,t') 



2tt poc 



drj / dt'T lc {0,u-^ m t - t'^P^^'-e' ,Lo'-t') 



2tt c og 



drj / dt'r ld (9,cj;(t),ir,t-t")P 1 ((t>,r),t";6',iv';t') 

OG J t' 



(78) 
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is the one loop contribution. The kernels ri Q , T\ c , and _ __§r U) — ui') n = ±1 (79) 

r 1( j are explicitly computed in Appendix lAl 2irN 

The expressions for the one loop contribution to T are 
rather complicated but several key features can be ex- 
tracted, namely 1) the introduction of a diffusion opera- 
tor, 2) a shift in the driving frequency, and 3) the addi- a± 2 (s;uj)P 1 (n,uj,uj' , s) 
tion of higher order harmonics to the coupling function K . . f ~ , 
/. The diffusion operator has the effect of shifting the ~ y#M J dv(3± 2 (s;uj,v)Pi(n,v,uj , s) 
marginal spectrum from the imaginary axis into the left i 
hand plane. The effect is that the finite size fluctuations = o ~ U U = 
to order 1/N stabilize the incoherent state. 

We can see these effects more easily by considering the 
special case of f(9) = sin# and g(u>) being a Lorentz dis- an< ^ 
tribution (see Appendix IB]). The Fourier-Laplace trans- 
formed equation of motion for the one loop propagator , \ fw / \ 1 r/ >\ i i r> 
i ,, n a n (s;id)Fi(n,LU,u> ,s) = - — —o(uj — u) ), \n > 2 
has the form nK ' 1 v ' ' ' ; y ' 11 



(81) 

g(uj)(l + /3i(s; w)) / dvPi(n, v, oj',s) where 
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(82) 



We can solve for Pi using the same kind of self- 
consistency computation that we used for Pq. This pro- 
duces an analogously defined dialectric function, A* (s). 

Ai(-) = 1 - f / (83 ) 

Stability of the incoherent state is determined by the 
spectrum of the operator T. Analogous to tree level, the 
continuous spectrum is given by the zeros of a n (s; ui) and 
the discrete spectrum by the zeros of A,* (s) given by (|83|) . 
However, at one loop order, the expressions for Pi rep- 
resent solutions to a coupled system of equations. Thus, 
the poles of tree level are shifted, and there are also new 
poles reflecting the interaction of the mean density with 
the two-point correlation function. These poles will have 
residues of 0(1/ N 2 ) owing to their higher order nature. 
For n = ±1,±2, we cannot solve for the poles exactly 



but we can approximate the shift in the tree level spec- 
trum by evaluating the loop correction at the value of the 
tree level pole, s = ^mw, which is equivalent to using 
the "on-shell" condition in field theory. Since the higher 
order modes will decay faster than the tree level modes, 
this essentially amounts to ignoring short time scales and 
is similar to the Bogoliubov approximation [25|, [26| . The 
remaining effective equation for Pi is now first order and, 
consequently, we can consider the spectrum of the implic- 
itly defined operator analogous to (|4ip . 

The continuous spectrum consists of all the zeros of the 
function a n . We expect a term of the form inui + 0(l/N) 
because of the tree level continuous spectrum. This will 
govern the behavior at large times. In this case, the "on- 
shell" condition is equivalent to Taylor expanding the 
loop correction via s — inu; + 0(1/ N) and keeping only 
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terms which are 0(l/N). This yields: 

a n (s; w) = a + in(u + Suj) + n 2 D 

where 



Suj 



K 2 



is a frequency shift and 



K > 



UJ' 



4 7 -K 
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D = 



A - 



2A 



,7-f) + 



(84) 



(85) 



(86) 



is a diffusion coefficient. The frequency shift, which is 
negative, serves to tighten the distribution around the 
average frequency. The diffusion operator serves to damp 
the modes which are marginal at tree level. 

The discrete spectrum arises from the zeroes of A* (s). 
We can again approximate the shift in the tree level zero 
by using the on-shell condition. This gives 



AiW = 1-f / 



duj- 



s + in(uj + Su>) 
A 2 2 7 - f 1 
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(87) 



We assume the shift will be small which allows us to write 
the zero of A* (s) as 
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7 ~T 



ds 



(88) 



where s = — (7 — -y) and <5A* (s) is the 0(1/N) correc- 
tion to A* (s). This results in 
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(89) 



Away from criticality (A <C 27) or for large A, this cor- 
rection is small. 

We conclude this section by writing down an effective 
equation of motion for the density function which incor- 
porates the effect of fluctuations. Recall the first equa- 
tion of the BBGKY hierarchy is 
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f(8' -8)p(x,t)p(x',t)d0'duj' 
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f(6' -9)C(x;x' ,t)d6'duj' '(90) 



This equation has a term (the "collision" integral) involv- 
ing the 2-point correlation function, C, on the right hand 
side. The equation determining the tree level propaga- 
tor (pn"j) is the linearization of the first BBGKY equation 
with C considered to be zero. The one loop correction 
to this equation ([75]) incorporates the effect of the cor- 
relations on the linearization. The diagrams in Figure [3] 
provide the linearization of the collision term, where C 
is considered as a functional of p. Using our one loop 
calculation, we can propose an effective density equation 
at one loop order 
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(91) 



where SI = uj + Suj and D is given by Eq. (f86]) . The field 
p{9, SI, t) is now defined in terms of the shifted frequency 
distribution 

G(Sl)^g(uj)(l-^-) (92) 

The new coupling constant K2(uj) is 0(1/N) and is due 
solely to the fluctuations. It arises from the term j3±2 
in the equation for P±. In fact, given the structure of 
the diagrams, it is clear that for 0(1 /N n ) there will be a 
new coupling, K n+ i, which corresponds to a sin[(n+ 1)9] 
term. In the language of field theory, all odd couplings 



are generated under renormalization. 

The generation of higher order couplings is especially 
interesting in light of the results of Crawford and Davies 
concerning the scaling of the density 77 beyond the onset 
of synchronization, i.e. r/ — po ~ (A — K c )@ [3l|, [Uj]. 
Although our calculations pertain to the incoherent state, 
the fact that the loop corrections generate higher order 
couplings is a general feature of the bulk theory defined 
by the action of Eq. (fT4|) as well. Thus, we expect a 
crossover from j3 = 1/2 to = 1 behavior to occur as 
N gets smaller. This is consistent with [3lj . wherein a 
crossover manifested as the rate constant became smaller 
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than the externally applied diffusion. In our case, the 
magnitude of the diffusion is governed by the distance to 
criticality and the number of oscillators. 

Our proposed effective equation (|9~Tj) is not self- 
consistent because we use the propagator to infer the 
form of the mean field equation. Thus, we neglect non- 
linear terms which may arise due to the loop corrections. 
In addition, our calculation applies specifically to per- 
turbations in the incoherent state. There are likely other 
terms we are neglecting for both of these reasons. The 
consistent approach would be to calculate the effective ac- 
tion to one loop order and derive the equation for p from 
that. This would involve essentially the same calculation 
we have performed here, but for arbitrary p(8,ui,t) (i.e. 
we would need to solve (|3Tj) for the propagator in the 
presence of an arbitrary mean) . 



measurement of the propagator per equation (|7Tj) . We 
perform simulations of TV oscillators with f(8) = sin 6*. 
We fix 2% of the oscillators at a specific angle (8q = 0) 
and driving frequency (unless N = 10, in which case we 
fix a single oscillator; the plots with TV = 10 have been 
rescaled to match the other data). The remaining oscilla- 
tors are initially uniformly distributed over angle 9 with 
driving frequencies drawn from a Lorentz distribution. 
We measure the real part of Z\(t). This measurement 
allows us to observe the behavior of the modes which are 
marginal at tree level. 



V. NUMERICAL SIMULATIONS 

We compare our analytical results to simulations of sin- 
gle oscillator perturbations, since this provides a direct 



Equation (|73| gives the behavior of 8Z\ (t) with a single 
oscillator fixed at 0$ and ojo at time t = 0. Recall that 
Z\ = in the incoherent state, so that we expect 8Z\ ss 
Z\. To tree level 
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In other words, the initially fixed oscillator has phase #0 ~ w o^! 110 information is lost as the time evolution progresses. 
Incorporating the one loop computation gives 
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where we have ignored a term of amplitude 0(1/N 2 ); we are only considering the contributions coming from the poles 
described in the previous section. With the one loop corrections taken into account, we see that Z\(t) relaxes back 
to zero as t — > 00. In the simulations, we compute the real part of Z\(t) with 8q = 0. This gives 



Re(Zi(t)) = 
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(95) 



The special case of u>o = and 80 = gives: 



Zi(t) = 
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-,-Dt 
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(96) 



The imaginary part vanishes so that Re(Zi(t)) = Z\{t). 

We first compare our estimate of the diffusion coef- 
ficient D given by ([SB")) with the simulations. We plot 
the measured decay constant D of Z\ compared to the 
theoretical estimate of (IM1) for the long time behavior in 
Figure O These data only include values of K = 0.3K C 
and K = 0.5i^ c (K c — 2j — 0.1). Higher values of K did 



not yield good fits due to the neglected contributions to 
Z\. These decay constants arc obtained via fitting the 
time evolution of Z\ to an exponential for t > 200s. In 
both cases they behave as 1/N for large N as predicted. 
There is a consistent discrepancy likely due to round- 
ing error after simulating for such a long period of time 
(~ 30000 time steps). This error appears as a small de- 
gree of noise which further damps the response, hence the 
decay constants appear slightly larger in Figure This 
effect can be seen in Figures [S] and [7] as well. For large 
times, the simulation data consistently fall slightly under 
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FIG. 5: The large time (> 200s) decay constants with zero 
driving frequency for the perturbed oscillator. Lines are the 
predictions given by Eq. I86|l . Solid line and circles represent 
K = 0.03. Dashed line and boxes represent K = 0.05. 



the analytic prediction. Similarly, the data is noisier at 
large times. 

Figures [6] and [7] show the evolution of Z\ (t) over time 
along with the analytical predictions and the tree level 
result for K = 0.3K C and K = 0.5K C respectively. For 
K = 0.3K C , the prediction works quite well, with perhaps 
the beginning of a systematic deviation appearing at N = 
10 and N = 50 (there is a slight initial overshoot followed 
by an undershoot at larger times). This same deviation 
is more pronounced for K = 0.5K C , although the data 
follow the prediction quite well nonetheless. 

Consistent with our expectations from the loop expan- 
sion, as we move closer to criticality, i.e. the onset of syn- 
chronization, the results for K = 0.7K C and K = 0.9K C 
do not fare as well. Figurc[S]dcmonstrates a marked devi- 
ation from the prediction. We have not shown analytical 
results for the lower values of N because the deviation 
is so severe. The same holds true for all the results for 
K = Q.9K C , so that we have just plotted the simula- 
tion data in Figure [9] The general trend of approaching 
the mean field result still holds. The primary feature to 
take from these plots is that the fluctuations increase the 
decay constant. The closer to criticality, the more impor- 
tant the fluctuations and the faster the decay, hence the 
systematic undershoot which grows as one nears critical- 
ity. The fastest relaxation to the incoherent state appears 
at high K for a given N and at low N for a given K, ei- 
ther limit results in increased effects from fluctuations. 
It would be necessary to carry the loop expansion to two 
or more loops in order to obtain good matches with these 
data. 

In Figures [TU] and QTJ we plot the time evolution of 
Z\ (i) given that the favored oscillator has a driving fre- 
quency of Ld = 0.05. Note first that Z\{t) approaches 
the tree level calculation as N — > oo. The amplitude 
of the oscillation also shows the same deviation as the 
Loo = data, namely that of a slight initial overshoot of 
the one loop prediction followed by an undershoot. In 



FIG. 6: Zi(t) vs. t for various values of N and K = 0.3K C . 
Each graph shows N = {10, 50, 100, 500, 1000}. Note that as 
N — > oo the curve approaches the tree level value. From top 
to bottom: Black line represents tree level. X's and violet line 
represent N = 1000. Triangles and blue line represent N = 
500. Diamonds and green line represent N = 100. Boxes and 
red line represent N = 50. Circles and purple line represent 
N = 10. 




Tree Level 

O N= 10 

□ N = 50 

O N = 1DD 

A N = 500 

x N=1000 



FIG. 7: Zi(t) vs. t for various values of N and K = 0.5A" C . 
Each graph shows N = {10, 50, 100, 500, 1000}. Note that as 
N — > oo the curve approaches the tree level value. Symbols 
as in Figure [6] 



addition to this, we can see an increasing frequency shift 
as TV — -> 0. The data, prediction, and mean field results 
eventually become out of phase. For intermediate values 
of N, one can see that the one loop correction follows 
this shift, while for N = 10, the mean field, data, and 
one loop results each have a different phase. 

In the case of n > 2 it is easier to write down a complete 
analytical solution for the time evolution. With too = 0, 
6*o = we have, 




Dn A t 



(97) 
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FIG. 8: Zi(t) vs. t for various values of N and K = 0.7K C . 
Each graph shows N = {10, 50, 100, 500, 1000}. Note that as 
N — » oo the curve approaches the tree level value. Symbols 
as in Figure [6] 

This is compared with a simulation result in Figure 1121 
We see the same general trends as the previous graphs. 
At large N, the simulation follows the prediction quite 
well. For small N, the simulation seems consistently 
higher than the one loop prediction with K = 0.3K C . 
For K = 0.5K C , the prediction is again sufficiently sin- 
gular that we have not plotted N = 10. The deviation is 
already apparent for N = 50. 

VI. DISCUSSION 

Using techniques from field theory, we have produced 
a theory which captures the fluctuations and correlations 
of the Kuramoto model of coupled oscillators. Although 
we have used the Kuramoto model as an example system, 
the methodology is readily extendible to other systems of 
coupled oscillators, even those which are not interacting 
via all-to-all couplings. Moreover, the methodology can 
be readily applied to any system which obeys a conti- 
nuity equation. We derive an action that describes the 
dynamics of the Kuramoto model. The path integral de- 
fined by this action constitutes an ensemble average over 
the configurations of the system, i.e. the phases and driv- 
ing frequencies of the oscillators. Because the dynamics 
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FIG. 9: Zi(t) vs. t for various values of N and K = 0.9K C . 
Each graph shows N = {10, 50, 100, 500, 1000}. Note that as 
N — * oo the curve approaches the tree level value. Symbols 
as in Figure [6] 



of the model are deterministic, this is equivalent to an 
ensemble average over initial phases and driving frequen- 
cies. Using the loop expansion, we can compute moments 
of the oscillator density function pcrturbatively with the 
inverse system size as an expansion parameter. However, 
it is important to point out that the loop expansion is 
equivalent to an expansion in 1/N only because of the 
all-to-all coupling. A local coupling will produce fluctu- 
ations which do not vanish in the thermodynamic limit. 

Our previous work in this direction developed a mo- 
ment hierarchy analogous to the BBGKY hierarchy in 
plasma physics. This paper fully encompasses that earlier 
work. The equations of motion for the multi-oscillator 
density functions derivable from the action are in fact the 
equations of that BBGKY hierarchy. The calculation in 
Ref. [Tt} is in the present context the tree level calculation 
of the 2-point correlation function, given by the Feyn- 
man graph in Figure [2^,. With the BBGKY hierarchy, 
the calculational approach involves arbitrary truncation 
at some order, with no a priori knowledge of how this 
approximation is related to the system size, N. Here we 
show that this approximation is entirely equivalent to the 
loop expansion approximation. Truncating the hierarchy 
at the nth moment is equivalent to truncating the loop 
expansion at the (n — l)th loop for the Zth moment. The 
one loop calculation is performed in the BBGKY context 
by considering the linear response in the presence of the 
2-point correlation function. This would produce a more 
roundabout manner of arriving at our one loop linear 
response. One should also compare our one loop calcula- 
tion with the Direct-Interaction- Approximation of fluid 
dynamics; the path integral approach in that context is 
the Martin-Siggia-Rose formalism [27[. Another possible 
equivalent means of approaching this problem is through 
the Ito Calculus, treating the density as a stochastic vari- 
able and developing a stochastic differential equation for 
it. 

An important aspect of our theory is that it is di- 
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FIG. 10: As the previous figures, but with K — 0.3K C and FIG. 11: As the previous figures, but with K — 0.5K C and 

ujo = 0.05. Solid line represents tree level. Dashed blue line u>o — 0.05. Symbols as in Figure [TD] 
represents the one loop calculation. Circles represent the sim- 
ulation data. 
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FIG. 12: 5Z a (t) vs. t for K = 0.3A' C (top) and K = 0.5iC c 
(bottom). Symbols as in Figure [6] 

rectly related to a Markov process derivable from the 
Kuramoto equation. One can employ the standard Doi- 
Pcliti method for deriving an action from a Markov pro- 
cess to arrive at the same theory. Although the Ku- 
ramoto model is deterministic, the probability distribu- 
tion evolves in a manner indistinguishable from a funda- 
mentally random process. The stochasticity of the effec- 
tive Markov process is due to the distribution of phases 
and driving frequencies. In other words, it is a state- 
ment about information available to us about the state 
of the system. The incoherent state is a state of high en- 
tropy. The single oscillator perturbation is one in which 
we have gained a small amount of information about the 
system and we ask a question concerning our knowledge 
about future states. In the mean field limit for the sin- 
gle oscillator perturbation, we always know where to find 
the perturbed oscillator given a prescription of its initial 
state. In the finite case, our ignorance of the positions 
and driving frequencies of the other oscillators makes a 
determination of its future location difficult. Eventually, 
we lose all ability to locate the perturbed oscillator as it 



interacts with the "heat bath" of the population. Fur- 
thermore, this result should be time reversal invariant. 
Just as we have no way of determining with accuracy 
where to find the oscillator in the future, likewise we 
have no means of determining where it has been at some 
time in the past. To prove this statement in the con- 
text of our theory would require an analysis of the "time 
reversed" theory, obtained essentially by switching the 
roles of <p and ip. The relevant propagator for this time 
reversed theory will be the solution of the linearization 
of the adjoint of the mean field equation. Accordingly, 
this adjoint theory will have loop corrections which will 
damp the time reversed propagator as well. 

It is important to point out that our formulation ac- 
counts for the local stability of the incoherent state 
p(0,u>,t) = g(u>)/2ir to linear perturbations along with 
demonstrating the order parameters Z n approach zero. 
In mean field theory, there is the possibility of quasiperi- 
odic oscillations so that Z n = 0, i.e. the modes de-phase, 
while the incoherent state is marginally stable and infor- 
mation of the initial state is retained. Our work shows 
that in a finite size system, this is does not happen; the 
incoherent state is linearly stable. 

We have considered exclusively the case of fluctua- 
tions about the incoherent state below the critical point. 
Above criticality, a fraction of the population synchro- 
nizes. In this case, to analyze the fluctuations one may 
need to employ a "low temperature" expansion in con- 
trast to our "high temperature" treatment. In essence, 
one separates the populations into locked and unlocked 
oscillators and derives a perturbation expansion from the 
locked action. At criticality, each term in the loop expan- 
sion diverges. This is an indication that fluctuations at 
all scales become relevant near the transition and thus 
a renormalization group approach is suggested. Our for- 
malism provides a natural basis for this approach. 

In summary, we have provided a method for deriving 
the statistics of theories defined via a Klimontovich, or 
continuity, equation for a number density. This method 
produces a consistent means for approximating arbitrary 
multi-point functions. In the case of all-to-all coupling, 
this approximation becomes a system size expansion. We 
have demonstrated further that the system size correc- 
tions are sufficient to render the incoherent state of the 
Kuramoto model stable to perturbations. 
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APPENDIX A: ONE LOOP CALCULATION OF THE PROPAGATOR 



The loop correction Ti a applied to P is given by 
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This term arises from one vertex with a single incoming line and two outgoing lines and one vertex with two incoming 
lines and a single outgoing line, hence the product of two tree level propagators, Pq. We can represent Ti a in 
Fourier/Laplace space as 
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If f(6) is odd then /(to) = —/(—to). Therefore. 
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In this form it is easy to see the different channels which appear in the correction. 
Evaluating this expression using the tree level propagator (|40[) gives us 
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The diagram Tic is given by 
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In Fourier space, we can write: 
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duj' 2 (imn)(m + n)f(—m — n) f (m) f (—m)P(n + m, uj' 3 , t; v, t')P(—m, uj, t; ui\, ti)P(m, ui' 2 , t'; uj[, ti) 
+ J duj' 2 inm(m + n)f(m)f{m)f(—m)P(—m 1 uj' 3 , t; u>\, ti)P(n + m, uj, t; v, t')P(m 1 cu' 2 , t'; oj[, t\) 
+ J duj 2 inm(m + n)f(—n — m)f(m)f{—n)P{n + m, u' 3l t; ui 2 , t')P(—m, U), t; LOi 1 t\)P(m 1 ui 2 , t'; ui[, ti) 
+ J duj 2 inm(m + n) f (to) f (m) / (— n)P(— m, uj' 3 , t; ui\, t\)P(n + to, u>, t; uj 2 , t')P(m, ui 2 , t'; uj[, ti) 
and taking the Laplace transform: 

t lc (n,uj,u,s) = (^-}2N 2 K 3 {2n) 3 J2 [ du' a du 1 du[daig(ui)g(w' 1 ) 

du)' 2 (imn)(m + n)f(—m — n)f{m)f(—m)P{n + m, uj' 3 ;v,s — s\)P(—m, u>;u>i, si)P{m)oJ 2 ;oj' 1 , — si) 
+ / du> 2 inm(m + n)f(rn)f(m)f(—rn)P(—m,u)' 3 ;u)i, si)P(n + m,u>; v,s — si)P(m, u)' 2 ; ui' x , — si) 
+ J duj 2 inm(m + n)f(—n — rn)f(m)f(—n)P(n + m, oj' 3 ; ui 2 , s — st)P(—m, uj; lo\, s\)P{m, u> 2 ', u)[, — si) 
+ J dcj 2 inm(m + n) / (m) / (m) / (—n)P(—m, u>' 3 ; u>\, Si)P(n + m, uj; uj 2 , s — si)P(m, lu 2 ; uj[, — Si) 
where the contour for si lies between and < s < — Re(s n ). Performing the integrals, we have: 
f lc (n,u),v,s) = T7 R ' 3 [J—j I dsi^2(imn)(m + n) 

f(-m - n)f(m)f(-m) 



1 



1 



A„_)_ m (s — si) s — si + i(rn + n)v s\ — imuj A_ m (si) 
1 



imKf(-m)/ VA m (-Si) 
1 



/(m)/(m)/(-m 
-1 



1 



imKf(m)J \A_ m (si) 
1 



imKf(—m) J \A m (— si) 
duj 2 f(-n - m)/(m)/(-n) 

-si + imuj 2 A TO (— si) 



1 



1 



1 (2TrN)P(n + m, uj; u, s - Si) 



A„+,„(s - si) s — si + i(m + ti)uj 2 s± - imuj A_ m (si) 



(A7) 



(A8) 



(A9) 
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+ / dw 2 f(m)f(m)f(-n) -. 



1 



imKf(m)J V.A_ m (si) 



- 1 (2irN)P(n + m, u; uj 2 , a - si) 



9(^2) 



1 



— si + imtj2 A TO (-Si) 



Finally, the diagram Tid is given by 

,.2/r 



00 />i 

/ dt'T ld (9,cj;<f>, v;t - t')P{<j>, v;9' ,u';t') 
■00 Jo 

2 



-N 2 



K 2 d 



d9' 3 duj' 3 dt 2 d9' i duj' i d9 i djjj l 



(2tt) 2 96> 
[/(^-%K). 9 K) 

{P (9' 3 ,u' 3 ,t;9 2 ,u 2 ,t 2 )P Q (9,uj,t;9 1 ,u; 1 ,t ) 
+ P (9' 3 ,u; 3 , t; 9 l ,u 1 ,t )Po{9, ^> *5 to, w a , t a )} 

^/(»2 ~ to) {Po(0 2 ,u' 2 ,t 2 ; 0[, u[, t )P(9 2 ,u 2 ,t 

+ P(9' 2 , J 2 , t 2 ; 9', <J, t')P (0 2 ,u 2 ,t 2 ; ^ , u[ ,t )}} 



Using the fact that / du'd9'P{9, w, t; 9', u 1 , f)g(w') = g{u)/N and / d9f{9) = Owc have 



(A10) 



(All) 



(A12) 



2tt pQO 



dp / dt?T u (9,u;<f>,v;t--e)P{<l>,v;6',w';lf) = 
00 Jo 

--^Y^q J d9' 3 dw' 3 dt 2 \\d9[dw' l d9 l dw l 

[f(9' 3 - 9)g(L> 1 )g(u>[)P (9 3 ,uJ 3 , t; 9 2 ,u 2 ,t 2 ) 
-^-f(9 , 2 -9 2 )P(9 2 ,u !2 ,t 2 ;9',u > ',t')] 



(A13) 



(A14) 



In Fourier-Laplace we have 



T ld (n;w,v;s) = —inKf(~n)g(uj) , 

N \A„(s) 



(A15) 



APPENDIX B: f(6) = sin (9 AND g(uj) LORENTZ 

In order to both simplify the correction and to provide 
a concrete example, we will specialize to the case that 
g(ui) is a Lorentz distribution and f(9) = sin9. f(9) = 
sin (9 is the traditional coupling for the Kuramoto model 
(and has the advantage of being bounded in Fourier space 
so that we avoid "ultraviolet" singularities) and using a 
Lorentz frequency distribution yields analytical results. 

The Lorentz frequency distribution is given by 



( \ 1 7 



7T 7^ + U> 



From this we can calculate 



A±i00 = 



s + 7 



K 



s + 7 



(Bl) 



(B2) 



and A„(s) = 1 for n ^ ±1. The residue of the function 
A±i(s) is 



Res 



A_ m (sx) 



(B3) 



and the pole is at s±i = — (7 — K/2). This provides a 
critical coupling 



A'c = 2 7 



(B4) 



above which the system begins to synchronize. The in- 
coherent state is reached when K < K c , which gives 
s± < 0. From f{9) — sin0 we also have /(±1) = 

In this case, diagram a evaluates to 
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f lo (n,w, v, s) 
A? 
N 



m=±l 



n(2m + n)f(—m) 



; — -—NP u (m + n,u)\ v, s + 7 ) 

2tt imKfim) 2 v ' ' ' ' 2 7 



n(2m + n)/(— m — n.) 
n(2m + n)/(— m — n) 



1 If 



1 



1 



s + 2 7 



2 



2-7T 2 - 7 - imz/) 27r ( s + 7 - -| + i(m + s + 27 - A' 

111 s + 7 — imv 

2ir k_ m (imv) 2ir (s — imv + i(m + n)uS) s + 7 — lp — im 



r 



(B5) 



(B6) 



There is no contribution for n = which we expect 
from probability conservation. The terms proportional 
to n(2m + n)f(—n — m) will always evaluate to 0, be- 
cause f(n ^ ±1) = 0. Since the tree level propagator 
contains such a term as well, we have the simplification 



T la (n,u>,v,s) 
N 4 

m=±l 



(B7) 



n(2m + n) 



#(a; — v) 



s + 7 — + i(m + n)w 



fi c (n,w,^, s) 



2iV 



For diagram c, we see that we can immediately ignore 
the third term because nmf(—n)f(—m — n) is always 0. 
We also see that the first term is only non-zero for n = 
±2, and the last term only for n = ±1. After performing 
the si integration, this leaves us with 



5(u) — v) 



1 



— > n(m + n) 

2^ V 1 s + 1 - K + i { n + m ) u 2 1 - A 

sM / 27-f \ 

+ inu \ 27 - A J 




2inuj 7 



K 



1 



|«W + 7 



s + 7 - 



-2»w + 7-f «+§»(" -a;) ( 7 _|) 
1 1 1 A 



s + 27 
n A 



A, 



7 - f + § 1/ 2 7 - A -7 + f - §iw 2 



2T 5M 2 7 -A s + 7 



s + 27 



f + f tw s + 7 



f + fii/ s + 2 7 - A 



(B8) 



Tirf is given simply by 



APPENDIX C: EQUIVALENT MARKOV 
PROCESS 



r ld (±l;a>,i/;s) 
firf(n 7^ ±1;cj,i/;s) 



= 



2 



1 A 



(B9) 



The action (fTBj) can be derived by applying the Doi- 
Peliti method to a Markov process equivalent to the Ku- 
ramoto dynamics. Consider a two dimensional lattice C, 
periodic in one dimension, with lattice constants ag in 
the periodic direction and a w in the other. (The radius 
of the eventual cylinder is R.) The indices i and j will 
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be used for the frequency and periodic domains, respec- 
tively. The oscillators obey an equation of the form: 



9i = v 



t,co 



(CI) 



The indices on 9 and lo run over the lattice points of the 
periodic and frequency variables. The state of the system 
is described by the number of oscillators m,j at each site. 
Given this, the fraction of oscillators found on the lattice 
sites is governed by the following Master equation: 



dP (n,t) 
~dt 



E 



13 

Vi,j-1 



-riijP(n,t) 

ag 

(mj-i + l)P(n j ,t) 



(C2) 



where the indices of the vector ft run over the lattice 
points. C, and fp (note the superscript) is equal to n 
except for the jth and j — 1st components. At those 
points we have n\ ■ = n^j — 1 and n 3 i j_ 1 = n^j-i + 1. 
The first term on the RHS represents the outward flux 
of oscillators from the state with rij oscillators at each 
periodic lattice point while the second term is the inward 
flux due to oscillators "hopping" from j — 1 to j. There 
is no flux in the other direction (u>); this lattice variable 
simply serves to label each oscillator by its fundamental 
frequency. 

Consider a generalization of the Kuramoto model of 
the form: 



(C3) 



N is the total number of oscillators and we impose /(0) 
0. The velocity in equation (|C1[) now has the form: 



N 



'* ,3 



(C4) 



i',3 1 



In the limit a 
The factor of n 



we have ia u = lo. Similarly, iag — > 9. 
v j has been added because the sum must 
cover all oscillators, and this factor describes the number 
at each site. We also sum over all frequency sites i' . The 
master equation (|C2[) now takes the form: 



dP (n,t) 



dt 



E 



xn. hj P{n,t) 



K 



"i ,3 



i',3 1 



K 



6 Na s ., 

x(n itj -x + l)P(fi j ,t)] 

= -HP{n,t) (C5) 
The matrix TL is the Hamiltonian. From this point, one 
can develop an operator representation as in Doi-Pcliti. 
Using coherent states and taking the continuum and ther- 
modynamic limits results in the action (fT"5|) . after "shift- 
ing" the field (p. 
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